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We consider the thermodynamically driven self-assembly of spheres onto the surface of a central 
sphere. This assembly process forms self-limiting, or terminal, anisotropic clusters (A/'-clusters) 
with well defined structures. We use Brownian dynamics to model the assembly of A/'-clusters 
varying in size from two to twelve outer spheres, and free energy calculations to predict the expected 
cluster sizes and shapes as a function of temperature and inner particle diameter. We show that 
the arrangements of outer spheres at finite temperatures are related to spherical codes, an ideal 
mathematical sequence of points corresponding to densest possible sphere packings. We demonstrate 
that temperature and the ratio of the diameters of the inner and outer spheres dictate cluster 
morphology and dynamics. We find that some A/'-clusters exhibit collective particle rearrangements, 
and these collective modes are unique to a given cluster size A^. We present a surprising result for 
the equilibrium structure of a 5-cluster, which prefers an asymmetric square pyramid arrangement 
over a more symmetric arrangement. Our results suggest a promising way to assemble anisotropic 
building blocks from constituent colloidal spheres. 



I. INTRODUCTION 

Anisotropic particles are compelling building blocks for 
self-assembled materials because their directional inter- 
actions can be exploited to create complicated and useful 
patterns [1-7]. One way to create anisotropic building 
blocks is to self-assemble them from simpler particles, 
where the building block represents a free-energy mini- 
mizing structure. Recently a number of papers have been 
published synthesizing and simulating compound build- 
ing blocks that are clusters of spheres[6, 8-18]. Colloidal 
spheres are attractive candidates for assembly because 
they can be made from a wide variety of polymers and 
metals, and their interaction potentials can be tuned with 
organic ligands, solvents, and salts. 

Here we consider a class of self- limiting, or "terminal" , 
colloidal clusters created by self- assembling a small pop- 
ulation of one type of particle, the "halo" particle (HP), 
around a second type of particle, the central particle 
(CP). The clusters are terminal because the only attrac- 
tive interaction is between the HP and CP, which are 
dilute in the fluid of HP, and therefore steric restrictions 
among co-adsorbed HPs inhibit further growth. The re- 
sulting clusters have structures determined by the in- 
teractions among the adsorbed HPs, which self-organize 
around the CP to minimize their free energy. 

Arrangements of HPs on the surface of a CP have been 
studied extensively by mathematicians in the context of 
optimal arrangements of points on a sphere [19-22]. The 
solutions provide a library of anisotropic clusters that 
can in principle be created with properly designed inter- 
actions among the constituent particles. In this work we 
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study hard sphere HPs that are attractive only to dilute 
CPs and not to other HPs, thereby producing clusters 
of HPs around a single CP. The arrangements of these 
HPs bear comparison to a particular set of solutions, the 
spherical codes, for certain ratios of particle diameters. 
We investigate the self-assembly of these clusters as a 
function of temperature, where entropy controls the equi- 
librium structure of the cluster, and in a semi-open sys- 
tem, where HP are free to bind and unbind from the CP 
surface. We also consider the effect of temperature on 
the cluster structures and dynamics at deviations from 
the perfectly dense packings that correspond to solutions 
of the spherical code. 

This paper is organized as follows. In Section II we 
briefly review sphere surface extremal point problems. In 
Section III, we introduce the methods we use to study the 
terminal A^-clusters, including Brownian dynamics sim- 
ulations, free energy calculations, and metrics for clus- 
ter structure and mobility. In Section IV, we report the 
results of our simulations, free energy calculations, and 
analyses. We find that terminal A/'-clusters self assemble 
across a range of diameters and temperatures and the 
structure of these clusters resemble spherical code solu- 
tions. These findings are supported by free energy cal- 
culations, which predict cluster sizes and distributions. 
Using Brownian dynamics and free energy calculations, 
we explain the surprising observation of a dominant low- 
symmetry N = b cluster, a deviation from the spherical 
code prediction. We calculate changes in cluster struc- 
ture across a range of diameter ratios and investigate the 
dynamics for different cluster sizes, including collective 
modes. We find that the dynamics for clusters of differ- 
ent sizes are different. In Section V, we discuss several 
ways this work can be extended to create more types 
of anisotropic particles via tuning of the particle inter- 
actions, constructing additional shells of particles, and 
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FIG. 1. (Color) A terminal A/'-cluster with an octahedral structure {N — 6) is self- assembled from a bath of HP and a CP. 
This cluster has applications as an anisotropic building block, could be used to manufacture a "patchy particle" by imparting 
patches on the CP at the contact points, or could be locked into a nanocolloidal cage structure. 



creating structurally reconfigurable particles. In Section 
VI we conclude with a summary of our findings. 



II. SPHERE SURFACE EXTREMAL POINTS 
AND SPHERICAL CODES 

The problem of finding extremal points obtained by 
optimally distributing points on the surface of a sphere 
to minimize a function / has been well studied in the field 
of mathematics [2 1-23]. The problem is typically posed 
as follows: 

Given N points on the surface of a sphere of radius R, 
what arrangement of the N points minimizes a function 

/? 

li f = k ^i^j r^J^ where rij is the Euclidean distance 
between the points i and j, and n = 1, minimizing / cor- 
responds to the Thomson problem, whose solution de- 
scribes the distribution of identical point charges on the 
surface of a sphere. As n ^ oo, the problem corresponds 
to the spherical code^ (also known as the Fejes Toth, or 
Tammes problem), whose solution maximizes the min- 
imum distance between any two sets of points [19-22]. 
Other possible choices for / include minimizing the max- 
imum distance of any point to its closest neighbor, also 
known as the sphere covering problem^ and maximizing 
the volume of the convex hull of the points. For each of 
these problems solutions are exactly known for some val- 
ues of while various numerical searches have suggested 
best solutions for other TV. For the functions mentioned, 
tables of putative solutions up to at least N = 130 can 
be found in Ref. [20]. 

Fig. 2 depicts the spherical code solutions for 1 < TV < 
12. The arrangement of points for TV = 4 corresponds 



to the vertices of a regular tetrahedron, TV = 6 an oc- 
tahedron, TV = 8 a square anti-prism, and TV = 12 an 
icosahedron. The point arrangement of TV = 11 is equal 
to the TV = 12 solution minus a single point, or an icosa- 
hedron with one truncated pentagonal face. For each TV, 
the point group - the group of isometrics that keeps one 
point fixed - of the arrangement [21] is shown in the up- 
per right corner. Each optimal arrangement of TV points 
on the surface of the sphere is unique except for TV = 5 
which has a continuum of solutions ranging from a trian- 
gular bipyramid (point group Dsh^ shown in Fig. 2 and 
Fig. 10b) to a square pyramid (point group C^y^ shown in 
Fig. 10a). All solutions in the continuum have two points 
at opposite poles of the central sphere and differ by the 
positions of the three remaining points on the equator. 
The square pyramid arrangement is equal to the TV = 6 
solution minus a single point. We discuss these structures 
in detail in Section IV D. 

If the TV points represent sphere centers, the spherical 
code solution corresponds to the densest packing of TV 
hard halo spheres that all "kiss" a central sphere. For 
any packing of spheres around a central sphere, we de- 
fine A to be the ratio of the central sphere diameter, Dc, 
to the halo sphere diameter, Dh- We denote the mini- 
mal possible diameter ratio for TV spheres, which corre- 
sponds to the spherical code solution, as Aat. In Fig. 2, 
An of each arrangement is shown to four significant dig- 
its in the bottom right corner. Notably, AAr=5 = AAr=6 
and AAr=ii = A7v=i2- In one of mathematics' most fa- 
mous debates, Isaac Newton and David Gregory argued 
whether the kissing number of unit spheres (A = 1) is 12 
or 13. Had it been known that a central unit sphere can 
only be kissed by 13 spheres if their radii is r < 0.9165, 
or Ai3 = 1.0911 [20], this would have settled the question. 
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tailed structural and dynamic quantities for each cluster. 



A. Hard Sphere and Sticky Sphere Model 

In a semi-open system, the spherical code solutions 
of Section II correspond to perfectly hard spheres ad- 
sorbed on a perfectly sticky sphere. Mathematically, 
perfectly hard spheres are points interacting via a func- 
tion that steps from infinity to zero (Fig. 3a) and per- 
fectly sticky spheres are points interacting via the same 
function plus an infinitely narrow square well function 
(Fig. 3b). In this work we use radially-shifted Weeks- 
Chandler-Andersen (WCA) and Morse models of hard 
and sticky spheres, respectively, which allow for compu- 
tational efficiency as well as direct comparison with their 
ideal mathematical counterparts (Fig 3). They also cap- 
ture, in a general sense, the repulsive and attractive in- 
teractions of the constituent particles we have in mind. 
As nanoparticle synthesis continues to mature, the types 
of interactions that can be used to guide the self-assembly 
of small particles can be precisely tuned over wide ranges 
of length and energy scales, and the models used in sim- 
ulations can be suitably adjusted. 
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FIG. 2. (Color) The arrangement of points (pink) that cor- 
respond to each spherical code solution for 1 < N < 12. The 
point group of each arrangement is shown to the upper right 
of each arrangement, and the densest packing diameter ratio 
Dc/Dh = An is shown to the lower right. For = 5, the 
triangular bipyramid configuration is shown. Other N = 5 
configurations are shown and discussed in Figures 10-9 . 
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Isaac Newton's conjecture that the kissing number is 12 
was not proven until 1953 [24]. 

We also note that the spherical code solutions for N = 
3-12, except = 5, are rigid or jammed. They contain 
no "rattlers", defined as spheres not in isostatic contact 
with other spheres [25, 26], and cannot be deformed other 
than global isometrics [25, 26]. 



III. METHODS 



(b) 3 




To predict and compare the terminal A/'-clusters of halo 
particles bonded to central particles we use computa- 
tional tools that sample equilibrium statistical mechan- 
ical ensembles. In particular, Brownian dynamics simu- 
lations of model particles are used to perform computer 
experiments of self-assembly and the results of these sim- 
ulations are compared against cluster probabilities calcu- 
lated from a free energy analysis based upon numerical 
partition function calculations [7]. We also calculate de- 



FIG. 3. (a) A mathematically ideal hard particle interaction is 
shown in solid black compared to the hard particle interaction 
(in dashed blue) given by the WCA potential (Eqn. 1). (b) A 
sticky sphere with a kissing contact potential when 5 ^ is 
compared to a model sticky sphere (in dashed blue) given by 
the Morse potential (Eqn. 3). 



The radially-shifted WCA potential is given by [27] 
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r > ^cutoff 



•(1) 



The shifting parameter a is defined as a = <j/i — cr, where 
CT/^ is the WCA "diameter" of the HP, and r cutoff = 
+ The interaction between two HPs can be made 
arbitrarily hard relative to their size by increasing ah- 
The cost of increasing ah is that the dimensionless time 
T that elapses over each time step is reduced as r oc 1/ah- 
We choose ah = 3cr for its computational efficiency and 
its relatively "hard" modeling of spheres. The energy 
parameter e also determines the "hardness" of the HP, 
as a larger energetic penalty to overlapping corresponds 
to a "harder" potential. The cost of increasing e is that a 
smaller simulation time step is needed to model a steeper 
function. The energy parameter is set to e = (0.1/T), 
where T is the temperature of the simulation, so that 
the hardness of the HP is independent of temperature. 

Using Eqn. 1 to model hard HPs is effective because 
of the large potential energy penalty associated with two 
HPs approaching closer than 3a. However, due to the 
soft nature of the potential, spheres with sufficient ki- 
netic energy can, in principle, approach as close as 2a. 
It is therefore useful to determine the effective hard par- 
ticle diameter of HP modeled by Eqn. 1. We use the 
Barker-Henderson equation [28] to calculate the effective 
diameter 



Dh,e 
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)dr 3.0786cr 
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where ^ = l/Zc^T, u{r) = UwcA and the potential 
is zero at asH = + a. For the purpose of as- 

sessing the error in our calculations based on the ef- 
fective diameter, we can characterize two HPs as con- 
tacting when the interaction energy between them is in 
the range < UwcA < lO/c^T, which corresponds to 
3.0a < Dh,e < 3.1225cr. 

The radially-shifted Morse potential [29] used to model 
the "kissing" contact potential between the HP and CP 
is given by 
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where Eq = b determines the depth of the energy well, 
P = 5.0/(7 determines the width of the energy well, and 
ro determines the radial displacement of the bottom of 
the energy well. The Morse potential interaction range 
is truncated at r cutoff = 2.5(7 + (ro — a). 

An effective CP diameter can be calculated by defining 
an HP and CP as bonded when the distance between the 
two particle centers is at the minimum of Eqn. 3. More 
properly, an HP and CP are bonded when they remain 
positionally correlated because the HP remains within 
a given displacement of the CP. We use the minima of 



Eqn. 3 and the effective diameter of Eqn. 2 to define an 
effective CP diameter Dc^e = 2ro — Dh^e- The ratio of the 
CP to HP diameter is therefore = Dc^e/Dh^ei foi" 
molecular dynamics). We choose to keep the hardness of 
the HP-HP interaction constant for all the MD simula- 
tions by holding Dh^e (i-e. ah) constant while varying ro 
to change I^c,e- 

The non-infinitesimal potential well width of Eqn. 3 
permits the bond between the CP and HPs to stretch 
a small amount while remaining bonded. At some A"^ 
ratios, this stretching, though small, may be enough to 
accommodate an additional HP bond to the CP. It is 
therefore useful to define the CP bond- stretched effec- 
tive diameter Di)s,c,e = "^Rhc — Dh,e^ where Rhc is the 
average center-to-center distance between a bonded HP 
and CP measured in a simulation. The bond-stretched 
diameter ratio is defined as A^^ = Di)s^c,e/Dh,e {^s for 
bond- stretched). A^^ is always greater than A^. When 
the cluster is loosely packed, the difference between the 
two measures converges to zero. 



B. Brownian Dynamics 

To model mixtures of halo particles and central parti- 
cles assembling in a thermal bath we perform Brownian 
dynamics (BD) simulations, implemented in HOOMD- 
blue[30]. The natural units of this system are: the ef- 
fective diameter of the HP, Dh^e = 3.0786<j; the mass 
of a HP, m; and the depth of the HP-CP energy well, 
Eq. The volume fraction, (^, is defined as the ratio of the 
total volume of the HPs and CPs to the simula tion bo x 
volume, the dimensionless time is t* = t/ {Dh^e\/^/ Eq) , 
and the dimensionless temperature is T* = ksT/EQ. We 
use periodic boundary conditions. Each particle is sub- 
jected to conservative, random, and drag forces, and its 
motion is governed by the Langevin equation discussed 
further in [31-33]. We use a value for the drag coefficient 
7 = 0.726 m/t*. The same drag coefficient is applied to 
HPs in both the free and bound state. The conserved 
forces between particles are per Eqns. 1 and 3 above. 



C. Free Energy Calculations 

The relative probability of finding a particular cluster 
of N HPs bound to a CP can be predicted using free en- 
ergy calculations detailed in references [7, 34, 35]. For a 
given A, the partition function is defined by the appro- 
priately weighted sum over all possible configurations of 
TV HPs bound to a CP for TV = 1 to oo. The contribu- 
tions of the distinguishable microstates to the partition 
function are calculated numerically. 

The partition function is calculated assuming ideal 
hard spheres and sticky spheres (Fig. 3). Given HPs of 
diameter Dh and a CP of diameter Dc^ the interaction 
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potential between ideal HPs is defined as 
Uh 



-H{r) = I 



oo r < Dh 
r>Dh 



(4) 



and the interaction potential between an ideal HP and 
an ideal CP is defined as 

r oc r<{D^^Dh)/2 
UH-c{r) = { -Eo r = {D, + Dh)/2 (5) 
[ r>{Dc^Dh)/2 

We define A-^ = Dc/Dh (/ for free energy calculation). 
As in section III A, to vary A-^, is held constant (set 
to Dh e from Eqn. 2) and Dc is changed. 

If Ajv=M < A/ < A^v^M+i? then configurations of M 
HPs bonded to the CP minimize the potential energy 
and configurations with more than M HPs have infinite 
potential energy (zero probability). Configurations with 
fewer than M HPs bonded to the CP increase the en- 
tropy of the cluster. When ksT ^ Eq^ the free energy 
can be minimized by clusters with fewer than M HPs, be- 
cause the entropy gained by the remaining HPs on the CP 
balances the increase in potential energy. In the grand 
canonical ensemble, at a fixed A-^, the probability of ob- 
serving a particular cluster s is given by the Boltzmann 
distribution: 



Q^^-0{Us-I^N) 

Z 



(6) 



where Z = exp(— /3(C/s — fiN)) is the partition 

function and Us—fJ^N = NEq. Without loss of generality, 
we treat /i =0. (A non-zero fi will only induce a uniform 
temperature shift in our final results.) 

In practice, calculating Z exactly is difficult, but by 
assuming that only a small number of clusters contribute 
to Z[7, 34, 35], the relative probabilities of these clusters 
can be determined. As in [7, 34, 35], the degeneracy Qs 
can be written as a product of three independent terms, 
the translational, Zt, rotational, Z^, and vibrational Zy 
partition functions. The translational partition function 
is approximately equal for all the clusters because they 
are all small compared to the accessible volume, and thus 
contributes equally to the fls of each cluster. 

To calculate the rotational and vibrational partition 
functions for an A^-cluster, we first assume an equilib- 
rium configuration defined by N HPs in a spherical code 
configuration at a radial displacement of {Dc + Dh)/2 
from the CP. The rotational partition function is then 
calculated as = c^-^? where is a temperature- 
dependent constant that is the same for all the clusters, 
/ is the determinant of the moment of inertia tensor, 
and K is the symmetry number of the spherical code con- 
figuration under rotation. Each sphere is given a unit 
mass. The vibrational partition function is proportional 
to the product of the vibrational freedom, or freedom 
to rattle, of each sphere in the cluster. The vibrational 
freedom of each HP can be measured as the fractional 
area of the surface of the CP it has access to, subject 



to the restrictions imposed by its neighboring spheres. 
We approximate the vibrational area available to a given 
HP in a particular configuration by using a Monte Carlo 
numerical approach whereby new positions for the HP 
are randomly generated and accepted if the HP does not 
overlap another HP. The accessible vibrational area is 
proportional to the total number of accepted positions 
that are part of a contiguous area that includes the HP's 
original position divided by the total number of random 
trials. If A-^ = A^v, when the diameter ratio matches the 
spherical code ratio, then most, if not all, of the spheres 
in an A/'-cluster are jammed and have no vibrational free- 
dom. 

The free energy calculation is approximate, as it does 
not consider the contribution of collective modes of HP 
motion to Z^ which, in certain systems can help stabi- 
lize one configuration over another [36]. As we show in 
Section IV B, each cluster has a small A range, A > A^v 
where collective modes are not present, and only local 
rattling is observed. Outside this range, we expect some 
error in the calculation of relative probabilities to ac- 
cumulate. The benefit of this free energy approxima- 
tion is demonstrated by both its favorable comparison 
to predictions made by BD simulations and by its abil- 
ity to rapidly predict the entire phase diagram. Apply- 
ing a more computationally intensive method to perform 
an exact free energy comparison would be an interesting 
topic for future study. 



D. Structure and mobility measures of a cluster 

The HPs in an A^-cluster for A = A^v are confined 
to a unique TV spherical code solution and cannot re- 
arrange or even rattle for N = 3 to 12, excepting N 
= 5. For A >> Atv, the HPs are free to randomly ar- 
range on the surface of the CP. We aim to understand the 
structure and dynamics of the A/'-cluster between these 
extremes. We perform BD simulations of pre- assembled 
clusters wherein the HPs are restricted to the surface 
of the CPs, a constraint imposed during the integration 
of the equations of motion. This allows the dynamics 
of HP rearrangement to be isolated from the dynamics 
of assembly and disassembly and prevents any stretch- 
ing of bonds from influencing the structures observed. 
Similar to section HI A the CP diameter is defined as 
Dc^e = 2ro — Dh^e^ where ro is the fixed distance between 
the CP and HP centers and Dh^e is the same as defined in 
Eqn. 2. The CP to HP diameter ratio in the constrained 
system is thus A^ = Dc^e/E>h,e (c for constrained) and 
varied by changing tq. Ac is initialized such that the HP 
can be sparsely randomly distributed on the surface (i.e. 
A^ >> Atv), and then slowly decreased over the course 
of a simulation to a target A^. 

The angular displacement between two HP bound to 
the same CP, or ^ = ZACB is defined by the centers 
of two HPs A and B and the CP C. To characterize 
the structure of the cluster, the distribution of angular 



6 



displacements between pairs of HPs, n{0) for = [0, tt] 
are measured for a fixed TV and over all HP pairs every 
10^ time steps during a simulation with 10^ total time 
steps. The value of n{0) for a given N and A*^ represents 
the likelihood of finding an HP at an angle 6 relative to 
a given HP, and J n{0)dO = N — 1. n{0) is analogous to 
a pair correlation function. 

To characterize the dynamics we calculate the time 
scale over which is no longer correlated with itself. We 
define the mobility parameter r from 

C{e{t),0{t^St))=e-^' (7) 

where C{6{t)^0{t + 6t)) is the normalized angular auto- 
correlation function and t is time. In this work, r has 
units of 1/10,000 time steps. The more mobile an HP 
is on the surface of the CP, the more rapidly its angu- 
lar displacement with respect to other HP decorrelates. 
When the rate of decay of angular correlations is zero, all 
the HP in the cluster are fully caged. We only calculate 
T for clusters that display more than one distinct peak 
in their n{0) distributions so that position swapping can 
be distinguished from local rattling. The lower bound on 
the T measurement is 1.5 • 10~^ because below this the 
HP position swaps occur too infrequently over a 10^ time 
step simulation for accurate values of r to be measured. 
We calculate r as a function of N and AA^ = A^ — A^v, 
or the difference between the diameter ratio of the TV- 
cluster and the N spherical code solution ratio. Because 
the HPs in the simulation are not perfectly hard and are 
constrained to the CP surface, it is possible for meaning- 
ful measurements to be made when AA^ < 0. 



E. The calculation of A 

In this paper, to elucidate different properties of N- 
clusters, several calculation methods are used, necessi- 
tating four different ways to determine A = Dc/D^^ the 
ratio of the central particle diameter, Dc^ to halo parti- 
cle diameter, D^. Each way was chosen to best represent 
the effective diameters of the HP and CP in the particu- 
lar method. These As are comparable to each other and 
to the spherical code solution ratios, A at of Fig. 2. We 
indicate the calculate type by the superscript x of A^, 
where x G {m, 65, c, /}, where the m, Brownian (molecu- 
lar) dynamics; 65, bond-stretched; c, constrained; and /, 
free energy calculation, as defined above. 



IV. RESULTS 

A. Self-assembly and free energy of A/^-clusters 

Using Brownian dynamics we simulate the self- 
assembly of clusters as a function A"^ and at two different 
temperatures to investigate the effect of thermal noise on 



the distribution of stable terminal A/'-clusters. We com- 
pare these results to the known spherical code solutions 
and to free energy calculations. 

Brownian dynamics simulations of self-assembly are 
initialized by placing 1000 CPs on a cubic lattice, spaced 
so as to behave as independent systems. The lattice is 
embedded in a bath of HPs at a total volume fraction 
of (j) = 0.24 — 0.27. The bath contains a minimum of 
four times as many HPs per CP as the maximum cluster 
size observed for that A"^. We perform a total of 760 
simulations of 20x10^ time steps, with time step size 
At* = 0.00363 at low (T* = 0.02) and high (T* = 0.1) 
temperatures. With this set of simulations, we calculate 
the cluster size distribution as a function of A^. 

At the low temperature we observe that the cluster 
sizes are highly monodisperse as a function of A"^. In 
Fig. 4, the mean cluster size assembled at the low tem- 
perature, T* = 0.02 for 0.01 < A"^ < 1 is shown as a 
black solid line. Grey shading indicates the range of clus- 
ter sizes observed at a particular A^. Over this range 
of A"^ clusters are uniform in size except when A"^ is 
near a value where there is a transition from one mean 
cluster size to another. At these transitions, we observe 
a narrowly distributed mixture of cluster sizes; e.g., at 
A^ = 0.71 for the T* = 0.02 curve, we find equal num- 
bers of clusters containing 8 or 9 HPs. 

In comparison to the low temperature data, the clus- 
ters at high temperature are both smaller on average, and 
have a broader distribution of sizes as a function of A"^. 
In Fig. 4, we show the distribution of clusters assembled 
at high temperature, T* = 0.1 (red). The region shaded 
pink represents the range of cluster sizes measured at a 
given A^ at T* = 0.1. At A^ = 0.71 for the T* = 0.1 
curve, we now observe clusters of 5, 6 7, and 8 HPs. We 
also observe that the N = b and N = 11 clusters are not 
stable at any A'^ at low temperature but are present in 
the broader distribution of clusters at high temperature. 

To test the stability of the self-assembled clusters at 
low temperature, we perform a simulation wherein the 
diameters of the CPs in large pre-assembled clusters are 
slowly decreased. A single system with A^ = 0.9489 
is equilibrated for 20x10^ time steps at T* = 0.02, at 
which time every CP is bonded to 12 HPs. Subsequently 
L^ce is decreased at a rate of 4.833 x 10~^a/At until A"^ 
= 0.0101. As discussed in reference [37], this decrease 
in the diameter is slow enough that the system remains 
quasi-static, i.e. the system is approximately in equilib- 
rium. For this system, as A'^ approaches a transition 
ratio, 1-3 HPs detach from a given CP and re-enter the 
bath, until only two HP are bonded to each CP. In effect, 
this quasi-statically decreased I^c,e simulation disassem- 
bles the clusters as a function of A"^. 

If bond-stretching is taken into account, we find that 
at a low temperature (T* = 0.02) the A^-clusters self- 
assemble at the A ratio predicted by the spherical code 
solutions. In tightly packed clusters, bond stretching 
makes A^^ > A^. In the bottom plot of Fig. 4, the lowest 
A^^ at which a cluster of size N is observed for the self- 
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FIG. 4. (Color) Top: The N clusters that self-assemble as a function of A*" and temperature is shown. The average N of the 
self-assembled cluster at T* = 0.02 is shown as a black line. The maximum and minimum N in the simulation is shaded grey. 
The average A'^ of the self-assembled cluster at T* = 0.1 is a red solid line. The maximum and minimum N in the simulation 
is shaded pink. Example clusters self-assembled in simulation are shown. Bottom: Accounting for bond-stretching and the 
effective diameter of the HP, the lowest ratio where a cluster of size A'^ observed in the quasi-statically decreasing simulation 
(blue triangles) and for the self-assembled simulations (black circles) are compared to the spherical code predictions (pink star). 
Error bars for the quasi-static simulation ratios are generated from the contact range of two HP. 



assembled (black circles) and quasi-statically decreased 
(blue triangles) simulation data is shown and compared 
to the spherical code A^r ratio (pink stars). Blue er- 
ror bars indicate the A'''' ranges from quasi-statically de- 
creased simulations, generated by assuming that the true 
diameter of an HP is the limits of the contact range de- 
fined in section III A. Good correspondence between the 
predicted and measured ratios is observed when bond- 
stretching and the appropriate eflFective diameters of the 
particles is accounted for. 

We calculate the free energies of all clusters from N = 2 
to A?^ = 12 over a temperature range of 0.02 < T* < 0.2 
and diameter ratio range of 0.05 < A-^ < 1.09. In the 
bottom left plot of Fig. 5, we report a "phase diagram" of 
the most probable cluster at each combination of T* and 
A-^. The plots in the bottom right and top left of Fig. 5 
show data from an in-page slice of the phase diagram at 
the low and high temperature, T* = 0.02 and T* = 0.1, 
and directly compare it to cluster distributions from the 
BD simulation data of Fig. 4. For example, the three 
clusters and probabilities depicted in the upper right of 



Fig. 5 are from a single high temperature BD simulation 
with A^ = 0.46. 

We see that the free energy calculations support the 
findings of the BD simulations. At high temperature and 
at a given A, both show a decrease in cluster size relative 
to the low temperature data, as well as a broadening in 
the distribution of cluster sizes. Discrepancies in peak 
height and shape between the two predictions in Fig. 5 
are likely due to the soft sphere approximation, not ac- 
counting for the change in the contact energy or effective 
diameter due to bond-stretching, and also to neglecting 
the collective vibrational modes in the free energy cal- 
culations. However, the free energy calculation shows 
that most of the features of the BD simulations at higher 
temperatures can be attributed to the offsetting of the 
increase in potential energy (i.e. fewer bonded HPs) by 
the commensurate increase in vibrational freedom of the 
remaining bonded HPs. 

Consistent with the BD simulation data, the free en- 
ergy calculation also predicts that N = 5 clusters are 
not stable at low (T* < 0.06) temperatures. Spherical 



8 




FIG. 5. (Color) The distributions of cluster sizes as a function of temperature and A as given by the free energy calculation 
and the BD simulations are compared. Bottom left corner: phase diagram of the free energy prediction of the most probable 
cluster size. Lower right and upper left corners: in-page slices of the probability of finding each cluster size Pn as predicted by 
the free energy calculation and BD simulation at the high and low temperature. Upper right corner: the three most common 
clusters found in the BD simulation at the high temperature and = 0.46. 



code solutions indicate that the densest = 5 clusters 
occur at the same A as the densest N = 6 cluster. Thus, 
at low temperature, when the free energy is dominated 
by the potential energy term, the = 6 cluster is al- 
ways stable over an = 5 cluster. However, the free 
energy calculation predicts that at T* > 0.06 there is a 
A range where an = 5 cluster is the most probable 
cluster. This A range is observable in the high temper- 
ature BD simulation data. The stabilization of the A^ 
= 5 cluster over the A" = 6 cluster at higher T* arises 
from the non-negligible contribution of the vibrational 
partition function, the only term in the partition func- 
tion that significantly differs between the two clusters. 
The A" = 11 cluster is similarly predicted to be unstable 
at low temperatures but stable over the A" = 12 clus- 



ter at a higher temperature. The free energy calculation 
also predicts an entropic stabilization of the N = 7 and 9 
clusters over the A" = 8 and 10 clusters, respectively, at 
higher T* and "triple points" , at which the probabilities 
of three clusters (e.g. 4, 5, and 6; or 7, 8, and 9) are 
equal. 



B. Structure of A^-clusters 

We next consider how the structure of each A^-cluster 
changes as A > A^v- Clusters that have a large range 
of A over which their structure is ordered and stable are 
desirable targets for synthesis. We investigate the struc- 
ture and dynamics of the clusters by modeling HPs con- 
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FIG. 6. (Color) The distribution of angular displacements n(0) for each cluster. The n{0) shows a structural fingerprint 
particular to each cluster. 
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FIG. 7. Cluster mobility as a function of the AA"" = A"" - Aa^. 
Note that for the = 6 and N = 12 clusters, the HPs do 
not become measurably mobile for AA^ >> 0. At the other 
extreme, A^ = 5 and A^ = 10 are mobile for AA"" < 0. Each 
data point is extracted from a linear fitting with a goodness 
of fit > 0.95. 



strained to the surface of a CP, as described in Section 
HID for different N and A'^. 

In simulation, we observe that a cluster of size N 
generally exhibits three different dynamics over different 
ranges of A^. In the first range, each HP remains locally 
caged. Each HP of the A^-cluster can be assigned to one 
point of the N spherical code solution of Fig. 2 and that 
mapping remains invariant under the dynamics of the 
cluster. Like an atom in a crystal lattice, each HP rattles 
about its point. In the second range, each HP can almost 
always be assigned to one point of the spherical code so- 
lution, however the mapping does not remain invariant 
under the dynamics of the cluster. The HPs sporadically 
rearrange but are still generally found rattling about the 
points of the spherical code solution. In the third range, 
the HPs cannot be assigned to points of the spherical 
code solution and move freely on the surface. Between 
the second and third range, we suspect that there is no 
distinct measurable boundary, but simply an increasing 
likeliness of a cluster being in "transitional" states. Be- 
low we show how the two measures introduced in section 
HID capture the signature features of these ranges. 

In Figure 6, the distribution of angular displacements, 
n(^), is shown for 2 < < 12 HPs constrained to the 
surface of a CP at T* = 0.02. For each TV, n(6>), is shown 
for three different A^, corresponding to A at, Aat+i and a 
midpoint between the two. These three distributions are 
shown in order of increasing A^ from left to right. Note 
that for A/" = 2, A7v=2 is zero, as it is always possible to 
add a second HP to a CP with one bound HP, regard- 



less of CP size. In this case, we arbitrarily choose the 
minimum A'^ = AAr+i/2. 

For each cluster we observe a unique n{0) structure 
fingerprint that softens as A^ increases. For A/" < 4, 
each HP has one equidistant ring of neighbors, result- 
ing in n{6) having a single peak that broadens as A^ 
increases. For A/" > 4, each n{0) has multiple peaks. For 
N > 4 except N = 5 and N = 10, the first peak at A at 
is narrow and not connected to other peaks, indicating 
HPs are locally caged at their spherical code points [38]. 
The width of a peak is proportional to the rattling of 
a HP within its local cage. As A^ increases, the peaks 
broaden and eventually become connected. This broad- 
ening and overlapping is associated with the degradation 
of the well-defined structure by increased rattling and 
sporadic rearranging. In no case did we find any evi- 
dence of new structures emerging. We note that for the 
clusters N = 5 and N = 10 the peaks are not distinct at 
the smallest A^ considered. For all if A^ >> Aat, then 
the HPs sample uniformly random arrangements on the 
CP surface and the n{0) distribution is a cosine function 
of ^, truncated to zero when is less than the angular 
diameter of the HP (e.g. in Fig. 6, A/" = 2, n{0) is a trun- 
cated cosine function for each A^). In Fig. 6 for A^ > 2, 
insofar as the distributions are far from converged to a 
cosine function, we observe structure derived from the 
underlying spherical code solution over the entire range 
of A considered 



C. Mobility of A^-clusters 

We next consider the dynamics of the HPs on the CP 
surface. As described in section HID, we can measure 
how rapidly the angular displacements of the HPs decor- 
relate at a given A^. We call this measure the mobility 
parameter, r. When r = 0, the HPs are in the first 
dynamical range; that is, each HP is fully caged and 
the angular displacement between any two HPs does not 
decorrelate. When r > 0, but small, the HPs are in 
the second dynamical range. In Figure 7, the r of dif- 
ferent clusters are calculated as a function of increasing 
A^ relative to the ratio at which the cluster is predicted 
to assemble, AA^ = A^ — An- We observe that the size 
of the first dynamical range varies widely among clus- 
ters. Noticeably, the N = 6 and N = 12 clusters are 
not measurably mobile until AA^ is large. Note that in 
Fig. 6, the midpoint n{0) data of both N = 6 and A^ = 
12 still have distinct separated peaks. In contrast, the 
N = 11 cluster becomes mobile at a much lower AA^ 
than A^ = 12 cluster, despite having nearly the same 
spherical code solution and AAr=ii = AAr=i2. Comparing 
the A^ = 11 and N = 12 distributions in Fig. 6, at A^ = 
0.9021 the two clusters have nearly identical n{0) distri- 
butions. At A^ = 0.9489, A^ = 11 is measurably mobile 
(r = 1.5 • 10~^)) but still has distinct peaks in n{0) that 
are only slightly softer than that of N = 12. By A^ = 
1.095 the n{0) peaks are noticeably softened for the A^ 
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= 11 cluster relative to that of the N = 12 cluster and 
the peaks are connected. At this ratio, the mobility of 
the = 11 cluster is r = 0.19 while the N = 12 cluster 
is just measurably mobile (r = 4 • 10~^). In comparison, 
the = 10 cluster is mobile even at the ratio at which 
it first self-assembles. The rapid increase of the mobility 
as a function of increasing in Fig. 7 is consistent with 
the connected peaks and the rapid softening of the peak 
structure for A^ = 10 in Fig. 6. 

The fact that HP mobilities for a particular cluster 
depend upon the cluster's structure is not surprising. 
However, it is not obvious that clusters of different sizes 
should have such variation in the widths of the first dy- 
namical range indicated in Fig. 7. There is little correla- 
tion, for example, between the mobility of the A" cluster 
and the range of A over which the A^ cluster is stable 
in Fig. 4. We find that the HPs in the A" = 5 and 
A" = 10 clusters are never fully caged, while the HPs 
for A" = 6 and A^ = 12 are fully caged for a large range 
of The mobilities for A" = 7, 8, 9, and 11 lie be- 
tween these extremes. Note that N = 6 and A" = 12 
clusters have highly symmetrical spherical code point ar- 
rangements with octahedral and icosahedral structures, 
respectively. Their HP centers define the vertices of Pla- 
tonic solids with equilateral triangle faces. For A" = 7, 8, 
9, 10, and 11, the convex polyhedra defined by the centers 
of the HP have pentagonal (A^ = 11), square (A^ = 8 and 
10) or nearly square (A" = 7, 9, and 10) faces. We hypoth- 
esize that these non-triangular "defects" in the spherical 
code solutions are responsible for the increased mobility 
of these clusters by providing locations where the bar- 
rier to rearrangement is low. However, for such small 
systems, the rearrangements of HPs in a mobile cluster 
is more appropriately viewed as a rearrangement of the 
entire cluster rather than a localized rearrangement. 

As discussed above, in the second dynamical range, the 
HPs in an A" cluster can almost always be mapped to the 
points of the A" spherical code solution, but that mapping 
does not remain invariant. We examine the A" = 5-12 
clusters at the A^ at which each cluster is first observed 
to be measurably mobile per Fig. 7 to understand how 
the HPs in a cluster rearrange. We find (Fig. 8) that the 
A' = 5-11 clusters each have a single unique (discounting 
reflections or rotations) rearrangement, which permutes 
the HPs over the spherical code points. The A" = 12 
cluster has two unique rearrangements. Short movies of 
these rearrangements can be found online in the supple- 
mental material. For the A^ = 5 and A' = 11 clusters, 
which have structures equivalent to the A" + 1 spherical 
code minus a single point, a rearrangement consists of a 
single HP "hopping" a gap to the available A" + 1 point. 
(The reason the A^ = 5 is found in this particular config- 
uration is discussed in the next section.) The clusters A" 
= 6, 10, and 12 exhibit a permutation whereby a ring of 
HPs rotate relative to the cluster in a manner resembling 
a twist of a Rubik's Cube™. The clusters A' = 7, 8, 9, 
and 12 exhibit a permutation whereby the cluster "buck- 
les" into a new permutation of the spherical code points. 
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FIG. 8. (Color) The rearrangements of clusters N — 5-12. N 
— 12 has two rearrangements. 



We find that the addition of the rearranging action for 
A" = 5-12 is sufficient to make each cluster ergodic. That 
is, every possible assignment of each HP to the spheri- 
cal code points can be explored by the cluster with no 
inaccessible microstates. This ergodicity is shown, using 
group theory, in the supplemental materials. 

For clusters A^ = 6, 7, 8, 9, 10, and 12, the rearranging 
action is a collective motion of particles in the cluster. 
Although this entropic contribution is not considered by 
the free energy calculation in Section IIIC, the free en- 
ergy calculations compare well with the BD simulations, 
demonstrating that local rattling is more important than 
collective modes for some ranges of A. 



D. Breaking the degeneracy for N = 6 

1. BD simulations 

The A^ = 5 spherical code has a continuum of solutions 
ranging from the vertices of a square pyramid to a trian- 
gular bipyramid. For dense A" = 5 clusters at non-zero 
temperature, we seek the relative likelihood of the clus- 
ter adopting particular configurations from the solution 
continuum. For this, we construct an order parameter 
that can distinguish between different configurations in 
our BD simulations. 

All A^ = 5 spherical code solutions have two points at 
opposite poles of the central sphere and differ by the po- 
sitions of the three remaining points on the equator. The 



order parameter is constructed by, first, dividing the five 
HPs into "pole" HPs and "equator" HPs. The neigh- 
bor distances, or distance between each HP and the four 
other HPs is measured. HPs that do not have one neigh- 
bor distance that > 1.2 times the distance of the other 
three neighbor distances are "equator" HPs. Second, the 
"equator" HP that is closest to other "equator" HPs or 
has the minimum summed neighbor distances is selected 
and and its center is labeled A. 

Finally, an angle measurement is constructed in the 
plane of the equator as follows. The centers of the pair 
of "pole" HPs are labeled Pi and P2. The points A, 
Pi, and P2 define a plane ^i. The centers of the two 
remaining HP are labeled Ei and E2. The line through 
El and E2 intersects 5*1 at Es and n is the normal vec- 
tor to 5*1. A plane ^2 orthogonal to ^i is constructed 
from the point Es^ A, and A -\- h. The coordinates are 
translated and rotated so that A and Es are both on 
the ^-axis of S2 and A has x-y coordinates (0, ro), where 
ro is the distance between the center of an HP and the 
CP. The origin corresponds to the center of the CP. The 
points El and E2 are projected to the plane 5*2 and the 
angles (< 7r/2) to the x-axis of S2 is measured. The order 
parameter x is defined as this angle, sampled twice per 
configuration. Each angle pair uniquely specifies a con- 
figuration in the solution continuum. A perfect square 
pyramid configuration corresponds to two measurements 
of X = and a perfect triangular bipyramid configuration 
corresponds to two measurements of x = 7r/6 0.524) 
radians. Fig. 9a illustrates how the order parameter was 
constructed, and shows a sampling of the HP positions 
in the 5*2 plane from a simulation at A^ = 0.4. The red 
circles correspond to the triangular bipyramid positions. 

We performed BD simulations of clusters of 5 HP at 
A^ = 0.4142 and 0.400 with T* = 0.02. Two histograms 
are shown of the sampled x at the two ratios, 0.4142 and 
0.4 in Fig. 9a and 9b, respectively. The figures show that 
the degenerate continuum of = 5 spherical code so- 
lutions is broken by the introduction of thermal noise. 
Surprisingly, we find that the square pyramid is the pre- 
ferred structure, even over the more symmetrical trian- 
gular bipyramid. As the cluster is packed tighter, an even 
stronger preference for the square pyramid configuration 
over other configurations emerges. 



2. Free Energy 

To understand the preference for the square pyramid 
configuration in the BD simulation we use a free energy 
calculation, which elucidates the role of entropy in break- 
ing the degeneracy. The more symmetrical triangular 
bipyramid configuration is used as the reference configu- 
ration. 

The square pyramid and triangular bipyramid HP clus- 
ters are shown in Fig. 10(a) and Fig. 10(b). In Fig. 10(c), 
using a free energy calculation, the probability of observ- 
ing the square pyramid relative to the triangular bipyra- 
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FIG. 9. (Color) (a) The order parameter x is constructed by 
measuring the angle of the particles on the equator. Scattered 
points from a simulation overlay an image of an SP configura- 
tion. Red circles indicate the sphere centers of a TBP config- 
uration. In (b) and (c) the distribution of x sampled in from 
a BD simulation is shown as a function of the diameter ratio 
A"" = 0.4142 and 0.4 respectively. 
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FIG. 10. (Color) (Color) (a) The square pyramid (SP) and 
(b) the triangular bipyramid (TBP) N = 5 spherical codes. 
The jammed and unjammed kissing spheres in each configu- 
ration are colored dark grey and pink, respectively. The path 
that the unjammed spheres can follow is traced on the central 
sphere. For (b) the central sphere is transparent so the full 
path around the equator can be seen. In the graph at the 
bottom, at the low temperature, T* = 0.02, the preference 
for the SP (black solid) over the TBP (red solid) is evident 
as the HP diameter approach the limiting packing diameter. 
This preference (black and red dashed lines) is even stronger 
at high temperature, T* =0.1. 



mid is shown at two temperatures as the HP diameter, 
Dh^ approaches the diameter of spheres corresponding to 
the densest possible packing, Dh N=5- For dense clusters, 
we observe the square pyramid is always the most likely 
configuration at nonzero temperature. We find this pref- 
erence is because the square pyramid has the most vibra- 
tional freedom. In Fig. 10(a) and Fig. 10(b), the locally 
unjammed HP in a cluster at A = An=^ are colored pink 
and the HP that are locally jammed are colored grey. 
In the degenerate continuum of = 5 spherical code 
solutions, only the square pyramid has only one locally 
jammed HP, and thus, the highest vibrational freedom. 

In Fig. 7, the = 5 cluster becomes decreasingly mo- 
bile as the cluster is packed tighter, i.e. on decreasing 
AA^. Unlike clusters for other values of A/", as r ^ 0, 



the HPs in the N = 5 cluster become locally caged for 
entropic, rather than energetic, reasons. 



V. DISCUSSION 

There are a number of ways the sticky sphere assem- 
bly method described above can be extended to create 
interesting new species of anisotropic particles. 

For example, we can now ponder a more general ques- 
tion. Given a desired arrangement of points, what 
HP-CP interactions and HP-HP interactions will re- 
sult in self-assembly of the arrangement? The analo- 
gous mathematical question was posed by L.L. Whyte 
in 1952, ''What spherical arrangements [of points] pos- 
sess extremal properties of any kind?^'[22] Ideally, we 
seek HP-HP interactions and HP-CP interactions that 
self-assemble repeatable and desirable patterns of HPs 
on the CP. 

Cohn and Kumar [39] show that all potential energy 
functions of distance that are completely monotonic, such 
as inverse power laws, share a subset of universally op- 
timal solution configurations. If the function is strictly 
completely monotonic, then the universally optimal so- 
lution is also unique. For points on the surface of a 
sphere, the only known universally optimal solutions are 
[39, 40] N = 1-4, 6, and 12; that is, a single point, an- 
tipodal points, points forming an equilateral triangle on 
the equator, and tetrahedral, octahedral, and icosahedral 
arrangement of points. For our purposes, this means 
certain desired point arrangements (e.g. a ring of 12 
points distributed around the equator of a sphere such 
as modeled in reference [41]) are likely to be inherently 
difficult to achieve from HP-HP interactions. Restrict- 
ing themselves to isotropic pair potentials and identical 
particles, Cohn and Kumar [42] constructed separate de- 
creasing convex potential energy functions that have cu- 
bic (N=8) and dodecahedral (N=20) configurations as 
their minimum. Thus references [39, 40, 42] imply that 
to assemble certain clusters, it will be necessary to use 
more complicated HPs with carefully constructed poten- 
tials, including non-completely monotonic or anisotropic 
interactions. 

Using an alternative approach, complex clusters may 
also be possible by simply adding stages to the assembly 
process. For example, if, after the terminal A/'-cluster of 
Fig. 1 is created, the bath of HPs is replaced by a bath 
of new HPs coated with the same complementary mate- 
rial as the CP, a second shell of spheres can be added 
to the first. The structure of this shell will also depend 
on the entropy and energy of the cluster at a given tem- 
perature. If the HPs in the second shell preferentially sit 
in the interstices of the first shell, the polyhedron they 
form will be the dual of the polyhedron of the first shell. 
This can make new types of point arrangements possi- 
ble. For example, the dual of the octahedron is the cube. 
A cubic arrangement of eight points on the surface of a 
sphere is not found as a minimum among most common 
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spherical surface functions [20]. A second shell of HPs 
that preferentially assemble the dual of the first shell of 
HPs may be a physically more viable method of assem- 
bling a cubic arrangement of spheres without requiring 
the elaborately constructed HP-HP interaction potential 
of Cohn and Kumar [42]. 

The results presented here may also be used to guide 
the synthesis of reconfigurable TV = 5 clusters. As shown 
in section IV D above, a small change in the packing frac- 
tion of the = 5 cluster introduces a significant change 
in the structure of the cluster. Thus, changing the effec- 
tive diameter of the central particle by a modest amount 
induces a switch between a relatively isotropic disordered 
cluster and an anisotropic square pyramidal cluster. 

There may also be a correspondence between other 
mathematical sequences of points distributed on a sphere 
and terminal cluster assembly problems. For example 
in reference [11, 12], clusters of cones and spheres were 
shown to form unique and precisely packed clusters aris- 
ing from free energy minimization subject to a convexity 
constraint. The authors find a packing sequence identical 
to that obtained from minimizing the second moment of 
the mass distribution of a cluster of particles constrained 
to a convex hull. We observe that the packing sequence 
produced in [11] also bears a strong resemblance to the 
distribution of points on the surface of a sphere that max- 
imizes the convex hull [20]. The authors further showed 
that sequence successfully describes the polyhedral struc- 
tures formed by colloidal spheres self- assembling on an 
evaporating droplet [13, 14]. That work serves another 
example of the correspondence between mathematical so- 
lutions of extremal points on the surface of a sphere and 
cluster structures obtainable in experiments. 

Another interesting variant to consider is a shaped cen- 
tral particle, as was done in reference [43] for the Thom- 
son problem. Considering the packing of HP around a 
shaped CP may lead to novel clusters and is a generally 
unexplored problem. 

VI. CONCLUSION 

In this paper we have demonstrated that hard 
and sticky spheres can self-assemble into terminal N- 
clusters with interesting and, in some cases, unexpected, 
anisotropics. These clusters have predictable preferred 
structures that depend on temperature and sphere diam- 
eter ratio. We find that some clusters exhibit collective 
particle rearrangements, and these collective modes are 



unique to a given cluster size. 

If assembled directly from a bath at low temperature, 
certain cluster sizes (e.g. N = 4, 6, 12) form robustly, 
while other clusters occur only over small ranges with 
relatively mobile structures (e.g. N = 7,9,10) and still 
others cannot be formed at all (e.g. A/" = 5,11). A "multi- 
step" process that assembles the clusters from a bath at 
a higher temperature, removes the bath, and lowers the 
temperature may enable these hard-to-form clusters to 
be formed robustly as well. It may even be possible to 
adjust the effective diameter of the HP or CP as a step 
in the assembly process. Our free energy calculations 
and Brownian (molecular) dynamics predictions of clus- 
ter structure provide a guide for designing such a process 
for optimal yield of a desired cluster size with a well- 
ordered structure. Clusters fabricated in this way may 
find use as building blocks for subsequent self-assembly, 
as templates for manufacturing precisely placed circular 
patches on the surface of a spherical particle, creating 
nanocoUoidal cages, or fabricating reconfigurable parti- 
cles. 



A. Acknowledgements 

We acknowledge Oleg Gang and Alexei Tkachenko for 
discussions of related problems. We acknowledge Daphne 
Klotsa for her helpful comments on the manuscript. CLP, 
MM and SCG were supported by the U.S. Department of 
Energy, Office of Basic Energy Sciences, Division of Ma- 
terials Sciences and Engineering, under award DE-FG02- 
02ER46000, the U.S. Department of Energy Computa- 
tional Science Graduate Fellowship. EJ received sup- 
port from the James S. McDonnell Foundation 21st Cen- 
tury Science Research Award/Studying Complex Sys- 
tems, grant no. 220020139 and from a National Defense 
Science and Engineering Graduate (NDSEG) Fellowship, 
32 CFR 168a. SCG is also supported by the DOD/DDRE 
under the National Security Science & Engineering Fac- 
ulty Fellowship award No. N00244-09-1-0062. Any opin- 
ions, findings, and conclusions or recommendations ex- 
pressed in this publication are those of the author (s) and 
do not necessarily reflect the views of the DOD/DDRE. 
This work was also supported by the Non-Equilibrium 
Energy Research Center (NERC), an Energy Frontier Re- 
search Center funded by the U.S. Department of Energy, 
Office of Science, Office of Basic Energy Sciences under 
Award Number DE-SC0000989. 



[1] S. C. Glotzer and M. J. Solomon, Nature Materials 6, 
557 (2007). 

[2] S. Sacanna and D. J. Pine, Current Opinion in Colloid 

& Interface Science (2011). 
[3] M. R. Jones, R. J. Macfarlane, B. Lee, J. Zhang, K. L. 

Young, A. J. Senesi, and C. A. Mirkin, Nature Materials 



9, 913 (2010). 

[4] S. N. Fejer, D. Chakrabarti, and D. J. Wales, ACS Nano 
4, 219 (2010). 

[5] S. N. Fejer, D. Chakrabarti, and D. J. Wales, Soft Matter 
7, 3553 (2011). 



15 



[6] A. J. Williamson, A. W. Wilber, J. R K. Doye, and 

A. A. Louis, Soft Matter 7, 3423 (2011). 
[7] E. Jankowski and S. Glotzer, Journal of Physical Chem- 
istry B (2011). 
[8] D. Wales and J. Doye, Journal of Physical Chemistry A 

101, 5111 (1997). 
[9] F. Sciortino, A. Giacometti, and G. Pastore, Phys. 
Chem. Chem. Phys. 12, 11869 (2010). 
[10] S. Mossa, F. Sciortino, P. Tartaglia, and E. Zaccarelli, 

Langmuir 20, 10756 (2004). 
[11] T. Chen, Z. Zhang, and S. C. Glotzer, Proc Natl Acad 

Sci 104, 717 (2007). 
[12] T. Chen, Z. Zhang, and S. C. Glotzer, Langmuir 23, 
6598 (2007). 

[13] V. N. Manoharan, M. T. Elsesser, and D. J. Pine, Science 

301, 483 (2003). 
[14] E. Lauga and M. P. Brenner, Phys. Rev. Lett. 93, 238301 

(2004). 

[15] Y.-S. Cho, G.-R. Yi, S.-H. Kim, D. J. Pine, and S.-M. 

Yang, Chemistry of Materials 17, 5006 (2005). 
[16] Y.-S. Cho, G.-R. Yi, J.-M. Lim, S.-H. Kim, V. N. 

Manoharan, D. J. Pine, and S.-M. Yang, Journal of the 

American Chemical Society 127, 15968 (2005). 
[17] Y.-S. Cho, G.-R. Yi, S.-H. Kim, M. T. Elsesser, D. R. 

Breed, and S.-M. Yang, Journal of Colloid and Interface 

Science 318, 124 (2008). 
[18] L. Hong, A. Cacciuto, E. Luijten, and S. Granick, Nano 

Letters 6, 2510 (2006). 
[19] F. Toth, Regular Figures (The Macmillan Company, 

1964). 

[20] N. J. A. Sloane, with the collaboration of R. H. Hardin, 
W. D. Smith and others. Tables of Spherical Codes, Cov- 
erings, Maximum Volume of Convex Hulls, and Min- 
imal Energy Arrangements published electronically at 
www . research . att . coin/~iij as/. 

[21] T. W. Melnyk, O. Knop, and W. R. Smith, Canadian 
Journal of Chemistry 55, 1745 (1977). 

[22] L. L. Whyte, The American Mathematical Monthly 59, 
pp. 606 (1952). 



[23] J. R. Edmundson, Acta Cryst. A48, 60 (1992). 

[24] T. Aste and D. Weaire, The Pursuit of Perfect Packing 

(Institute of Physics Pubhshing, 2000). 
[25] H. Cohn, Y. Jiao, A. Kumar, and S. Torquato, ArXiv 

e-prints (2011), arXiv:1102.5060 [math.MG]. 
[26] T. Tarnai and G. Zs., Mathematical Proceedings of the 

Cambridge Philosphical Society 93, 191 (1983). 
[27] J. D. Weeks, D. Chandler, and H. C. Andersen, The 

Journal of Chemical Physics 54, 5237 (1971). 
[28] J. A. Barker and D. Henderson, Journal of Chemical 

Physics 47, 4714 (1967). 
[29] P. M. Morse, Phys. Rev. 34, 57 (1929). 
[30] "HOOMD-blue," http://codeblue.umich.edu/ 

hoomd-blue/ (2010). 
[31] Z. Zhang, M. Horsch, M. Lamm, and S. Glotzer, Nano 

Letters 3, 1341 (2003). 
[32] C. lacovella, M. Horsch, Z. Zhang, and S. Glotzer, Lang- 
muir 21, 9488 (2005). 
[33] C. lacovella and S. Glotzer, Journal of Chemical Physics 

129 (2008). 

[34] D. J. McGinty, J. Chem. Phys. 55, 580 (1971). 

[35] G. Meng, N. Arkus, M. P. Brenner, and V. N. Manoha- 
ran, Science 327, 560 (2010). 

[36] A. Haji-Akbari, M. Engel, and S. C. Glotzer, The Jour- 
nal of Chemical Physics 135, 194101 (2011). 

[37] C. L. Phillips and S. C. Glotzer, submitted. 

[38] E. R. Weeks and D. A. Weitz, Chemical Physics 284, 361 
(2002). 

[39] H. Cohn and A. Kumar, J. Amer. Math Soc. 20, 99 
(2007). 

[40] B. Ballinger, G. Blekherman, H. Cohn, N. Giansiracusa, 
E. Kelly, and A. Schiirmann, Experimental Mathematics 
18, 257 (2009). 

[41] Z. Zhang and S. C. Glotzer, Nano Letters 4, 1407 (2004). 
[42] H. Cohn and A. Kumar, ArXiv e-prints (2009), 

arXiv:0906.3550 [cond-mat.stat-mech] . 
[43] G. Vernizzi and M. Olvera de la Cruz, Proceedings of the 

National Academy of Sciences 104, 18382 (2007). 



